# Alexander F. Gazmararian
# afg2@princeton.edu
#January 9, 2024

#Load packages
library(tidyverse)
library(tidylog)
library(here)

# Load data
g <- readRDS(here("data", "inter", "acs_raw.rds"))
g$GEOID1 <- as.numeric(g$GEOID1)
g$GEOID2 <- as.numeric(g$GEOID2)

# Load matched counties data
ct <- readRDS(here("data", "output", "matched_df.rds"))
ct <- subset(ct, select = c(fips, treat))
ct <- unique(ct)

# Merge migration flows with treatment status
g <- left_join(g, ct, by = c("GEOID1" = "fips"))
g <- left_join(g, ct, by = c("GEOID2" = "fips"))

# Rename treatment variables
g <- g %>%
  rename(
    treat1 = treat.x,
    treat2 = treat.y
  )
names(g) <- tolower(names(g))

# Estimate flows from treated to other counties
g.all <- subset(g, age_label == "All ages")

g.all %>%
  group_by(geoid1, year, treat1) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  group_by(year, treat1) %>%
  summarise(movednet = mean(movednet, na.rm = TRUE))


g.all %>%
  filter(treat1 == 1) %>%
  group_by(geoid1, year) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  group_by(year) %>%
  summarise(movednet = median(movednet, na.rm = TRUE))
g.all %>%
  filter(treat1 == 1) %>%
  group_by(geoid1, year) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  group_by(year) %>%
  summarise(movednet = mean(movednet, na.rm = TRUE))

## Figure I2----
p.outflow <- g.all %>%
  filter(treat1 == 1) %>%
  group_by(geoid1, year) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  ggplot() +
  geom_histogram(aes(x = movednet)) + 
  facet_wrap(~ year) +
  theme_bw(base_size = 13) +
  labs(x = "Net Migration", y = "Count") +
  theme(panel.grid = element_blank()) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 50))
p.outflow

ggsave(
  p.outflow,
  filename = here("output", "figures", "si_fig_I2_netmigrate_all.png"),
  width = 6.5,
  height = 2.9,
  scale = 1.5,
  dpi = 300
)


## Figure I3----
# Estimate flows from treated to control counties
g.all %>%
  filter(treat1 == 1 & treat2 == 0) %>%
  group_by(geoid1, year) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  group_by(year) %>%
  summarise(movednet = mean(movednet, na.rm = TRUE))
g.all %>%
  filter(treat1 == 1 & treat2 == 0) %>%
  group_by(geoid1, year) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  group_by(year) %>%
  summarise(movednet = median(movednet, na.rm = TRUE))


p.flowtocontrol <-
  g.all %>%
  filter(treat1 == 1 & treat2 == 0) %>%
  group_by(geoid1, year) %>%
  summarise(movednet = sum(movednet, na.rm = TRUE)) %>%
  ggplot() +
  geom_histogram(aes(x = movednet)) +
  facet_wrap(~ year) +
  theme_bw(base_size = 13) +
  labs(x = "Net Migration", y = "Count") +
  theme(panel.grid = element_blank()) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 100))
p.flowtocontrol
ggsave(
  p.flowtocontrol,
  filename = here("output", "figures", "si_fig_I3_netmigrate_2control.png"),
  width = 6.5,
  height = 2.9,
  scale = 1.5,
  dpi = 300
)
